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Results of a numerical solution for radiative heat transfer in homogeneous and nonhomogeneous participating 
media are presented. The geometry of interest is a finite axisymmetric cylindrical enclosure. The integral 
formulation for radiative transport is solved by the YIX method. A three-dimensional solution scheme is applied 
to two-dimensional axisymmetric geometry to simplify kernel calculations and to avoid difficulties associated 
with treating boundary conditions. As part of the etTort to improve modeling capabilities for turbulent jet 
diffusion flames, predicted distributions for flame temperature and soot volume fraction are used to calculate 
radiative heat transfer from soot particles in such flames, it is shown that the nonhomogeneity of radiative 
property has very significant effects. The peak value of the divergence of radiative heat flux could be under- 
estimated by a factor of 7 if a mean homogeneous radiative property is used. Since recent studies have shown 
that scattering by soot agglomerates is significant in flames, the effect of magnitude of scattering is also inves- 
tigated and found to be nonnegligible. 
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Nomenclature 

= area, m : 

= absorption coefficient, 

= coefficient of scattering phase function, Eq. (2) 
= particle diameter, pim 
= emissive power of the medium, W/m- 
= radiation intensity, W/m : sr 
= kernel of integral equation. Eq. (8) 

= distance between boundaries, m 
= number of angular quadrature ordinates 
= inward unit normal vector 
= Legendre function of order n 
- radiative heat flux vector. W/m : 

= radiation quantity 

= distance from a point within the medium to the 
wall in the direction co t , Eq. (11) 

= position vector 
= cylinder radius, m 
= expansion function of Legendre series 
= scattering coefficient, m 
= temperature, K 
= volume, m 3 

= weights of discrete ordinates set 
= radiative heat flux, W/m : 

= distance, m 

= cylinder height, m . . 

= coefficient of anisotropic scattering phase 
function 

= surface emissivity 
= direction cosine 
= extinction coefficient, m" 1 
= direction cosine 
= direction cosine 

= density of the gaseous mixture, kg/m J 
= Stefan-Boltzmann constant. 

5.6696 x 10' H W/m : K J 
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= optical thickness, xrr () 

= normalized scattering phase function 
= domain or boundary of integration, Eqs. (5-7) 
= scattering albedo; solid angle 
= unit vector, Eq. (9) 


Sub scrip is 

b ~ blackbody 

g - participating medium for radiative transfer 

r = radiation 

s — boundary surface 


Superscript 

f = incident direction 


Introduction 

R ADIATIVE heat transfer in a cylindrical enclosure with 
a participating medium is a problem of practical impor- 
tance, e.g., in the design of industrial furnaces and many 
combustion devices. A solution method that is accurate, ef- 
ficient of both computing time and storage, flexible in a com- 
plex geometry, and compatible with the energy equation is 
needed for the prediction of radiative performance. Few cur- 
rent methods can satisfy all or part of these requirements. 
The Monte Carlo method is flexible and requires little storage, 
but can be extremely time-consuming, and the results are 
subject to statistical error. The zonal method and the finite 
element method are accurate but require large amounts of 
computer time and storage. It is also difficult to handle non- 
homogeneity or anisotropic scattering using the zonal method. 
The product integration method, 1 while faster than the zonal 
and finite element methods, does not reduce the required 
storage. The discrete ordinates Sn method, although accurate 
and less demanding of memory for large grid systems, suffers 
from ray effects 2 and high computer time for multidimensional 
combined mode heat transfer problems. To deal with multi- 
dimensional nonhomogeneous media, adaptive grid and adap- 
tive difference schemes must be used with the discrete-ordi- 
nates method to maintain the same order of local error. This 
will consume large amounts of CPU time and memory, and 
make calculation using the (Sn) method impractical. The 
spherical harmonics (Pn) method, 3 which needs high-order 
approximation to achieve accurate results for an optically thin 
medium, is tedious in formulation and also requires large 
amounts of computer time and memory in dealing with non- 
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homogeneous media. Recently, the finite volume method has 
been applied to cylindrical geometry. 4 For problems involving 
fluid flow computation, the method can be formulated to use 
the same grid system. However, a rather involved discreti- 
zation scheme must be used. 

In this article, a recently developed numerical method, the 
YIX method 5 - 6 and its extension to three-dimensional ge- 
ometry, 7 is used to solve the axisymmetric radiative transfer 
problem within turbulent jet diffusion flames. The formula- 
tion and solution scheme with cylindrical geometry (not lim- 
ited to axisymmetric) is described in detail. Results are pre- 
sented for homogeneous and nonhomogeneous participating 
media. The significance of the nonhomogeneity effect is dem- 
onstrated through comparisons with results for equivalent ho- 
mogeneous media. The treatment of media with spectrally 
dependent radiative properties is not considered in this study, 
but will be a subject for a future paper. 

There is increasing interest in the calculation of radiation 
heat transfer within nonhomogeneous sooting flames.*" 10 As 
part of the effort to couple turbulent diffusion flame modeling, 
soot formation and oxidation, and radiation heat transfer in 
a flame code, the treatment of radiation heat transfer within 
nonhomogeneous, absorbing, emitting, and scattering media 
is discussed. 


Mathematical Formulation 

The radiative heat transfer equation is written as 11 


d/(r, oj) 

~d/ 


- x/(r, oj) + ai b (r) 
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/(r, a/) dtr/ 


(1) 


where the phase function, based on the Mie theory, can be 
expressed as‘ : 


,v 

a/) = T ( 2 /t + \)a„P„(o>' w') (2) 


For particles with the size parameter (ttO/A) much less than 
unity, the scattering effect is usually neglected, e.g., Rayleigh 
particles. P,.(oj -w') can be expanded using the addition theo- 
rem of the Legendre function, 13 and 

u)-oj' = /x/x' + ff' + w' (2) 

where (/x, £, tj) are the direction cosines of oj with respect to 
each coordinate axis. After the terms are expanded, 

the phase function can be expressed as 

N 

a /) — ^ (2 n -I- \)a„P„{w’ a/) 

= 2 ( 4 ) 


where ft, = 1, /3, = /3 : = /3, = 3a, , ftt - 5a 2 /4, = /3 ft = 

5aJ3, 0 7 = ft< = 5a : /12, etc. 

The integral formulation of radiative heat transfer in a gen- 
eral three-dimensional, gray, emitting, absorbing, and ani- 
sotropic scattering medium corresponding to Eq. (1) by Tan 1 
is used here. Crosbie and Farrell 14 also developed similar 
integral expressions for intensity in three-dimensional cylin- 
drical geometry. The present formulation is convenient to 
couple with the energy equation, since heat flux and its di- 


vergence are computed directly in addition to the computa- 
tional efficiency for high-order scattering phase function 
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In the above equations, e K and e v are the blackbody emissive 
powers of the medium and the boundary, and q s is the net 
radiative heat flux at the wall. For isotropic scattering and 
nonscattering media, the second terms on the right sides of 
Eqs. (5-7) can be deleted. The kernel K is 


exp 


K(r, r') = 


[-r< 


r + ojt) dr 


( 8 ) 


and the unit vector 


When the medium is nongray, the integral equations are 
essentially the same, except that all radiative quantities are 
wavelength-dependent, and e g and e t are replaced by the spec- 
tral Planck function. 


Numerical Method 


Geometry 

Unlike planar two-dimensional (x-y) geometry, which ex- 
tends infinitely in the third ( z ) direction, axisymmetric cylin- 
drical (r-z) geometry does not extend infinitely in the third 
(0) dimension. For the planar two-dimensional problem, by 
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Fig. 1 Geometry of the cylindrical enclosure and coordinate system. 
Note the example ray of ab originated frorn the volume node a and 
hit boundary at point b. The direction of ab is determined from the 
discrete ordinates set used in angular quadrature. 

integrating with respect to the infinite length direction in the 
volume integration in Eqs. (5-7), the Bickley functions of 
different degrees will be the kernels. However, if one per- 
forms the volume integrations with respect to 6, as in the case 
of axisymmetry geometry (Fig. 1) 

I f f n K(r, r)w k (r)SM dV(r') 

- (* [ [ K(r, r')w k (r )S k (oo)r' d r d0' dz ' 

Jo Jn J ii 

depending on S k , the integration related to 0can be expressed 
as 


f 


exp 



x(r -h cot) dt 


d<9' 


a flat wall. 15 These difficulties are eliminated in three-dimen- 
sional schemes. 

Numerical Quadratures 

The integrations of Eqs. (5-7) are performed using the 
YIX method. 5 The integral equations are first rewritten into 
the distance-angular form. To maintain the same order of 
accuracy in angular integration at volume and boundary ele- 
ments in three-dimensional geometries, the fully symmetric 
discrete ordinates and* weights sets were used. 17 The use of 
discrete ordinates sets is discussed in Hsu et al. 7 

The volume and surface integrations on the right sides of 
Eqs. (5-7) are constructed as follows: 
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where the R(r , <*>) is defined as 
R(r , io) = min | r 4- cot | 

r - ojf./X) 

= length of a beam emitted from r in w direction 
and striking the nearest boundary 


p* expf - xVr : 4- /*'- — 2rr cos 6' + (z — z') : ] 

= — - - — —dd' 

Ji» ir[r 4- r’~ - 2 rr' cos 6' + (z — z') : ]"' 2 

(10) 

where r and f are any two points within the cylinder, and 0 ' 
is the angle between r and r projected on the x-y plane. 
Equation (10) is a kernel function different from the Bickley 
function. The kernel function is a function of three variables 
( r . r'. and z-z'), which makes it difficult to tabulate or ap- 
proximate by any computationally efficient function. Al- 
though direct numerical integration is possible, the amount 
of computation involved makes it essentially equivalent to a 
three-dimensional calculation. Another approach is to assume 
that the grid points and angular quadrature ordinates are 
predetermined, so that the new kernel function can be re- 
duced to a single variable. The YIX integration points can 
then be precalculated and stored. However, the significantly 
larger number of YIX integration points will also make the 
computation equivalent to a three-dimensional calculation. 

Other difficulties in treating axisymmetric cylinder cases as 
a two-dimensional problem are 1) an artificial symmetry 
boundary condition must be imposed at r ~ 0 (Menguc and 
Viskanta. 15 and Jamaluddin and Smith 1 *), which could result 
in singularity caused by 1/r-term 15 ; and 2) how to account for 
the curvature of the wall at r ~ r i)y rather than treating it as 


N w is the number of ordinates, which depends on the order 
of discrete ordinate sets used. For the Sn discrete ordinate 
set, N w = n(n + 2). The distance integrals in Eq. (11) are 
evaluated using the YIX quadrature. It is interesting to note 
that the right sides of Eqs. (5) and (6) are essentially the same 
except for the additional 5, term in Eq. (6). Therefore, the 
integrations in Eq. (6) can be avoided and a significant re- 
duction of computational time can be achieved. For calcu- 
lation with a high-degree anisotropic scattering phase func- 
tion, this causes only a minor increase in CPU time [due to 
the second term on the right side of Eq. (5)] compared with 
other methods. This is impossible for the two-dimensional 
formulation, where different kernels exist in Eqs. (5) and (6). 
Additionally, the current scheme is flexibile enough to use 
multiple discrete ordinates sets in the formulation, which is 
very advantageous in dealing with ray effects. 7 This flexibility 
does not exist in the conventional discrete ordinates method 2 -'* 
for the differential-integral formulation of the radiative trans- 
fer equation, where a consistent discrete ordinate set must be 
used even in the homogeneous case. By utilizing the axisym- 
metric geometry, only the first “wedge” of the boundary and 
volume elements is calculated (Fig. 1). To reduce computa- 
tional time further, the angular quadrature is not carried out 
for all N w rays: in practice, only the rays with positive y 
ordinates are calculated, due to the symmetry. The current 
scheme thus has the computational advantages of the three- 
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dimensional scheme, but at a cost approaching that of the 
two-dimensional scheme. 

Solution Procedure 

The discretized integrals of Eqs. (5-7) are solved by iter- 
ation. The steps are 1) give an initial guess for V *q, % iv*(r), 
and qX r )* with e A r )' an ^ e si r ) known; 2) calculate integrals 
on the right side of Eqs. (5-7) by the YIX quadrature; 3) 
calculate the new V q rf w t (r), and q x (r)\ and 4) go to step 2 
unless the convergence criterion is satisfied. 

For the case with e k .(r), and e t (r) as the unknowns, V q/\s 
obtained from the given heat source term of the energy equa- 
tion and q s (r) from the given flux boundary conditions. The 
iterative procedure is the same as above. 

Results 

The computation was performed on a Sun 4/690 worksta- 
tion. In all calculations, the first YIX integration point is 0.01 , 
except in the optically thin case, with overall optical thickness 
less than 0.1, where a first integration point of 0.001 is used. 
Typical run times can be several seconds to several hundred 
seconds, depending on the solution accuracy required and the 
optical thickness of the problem. Run time will be discussed 
later for the examples shown in Figs. 2 and 3. 

In Fig. 2, the YIX method is applied to a benchmark prob- 
lem: a uniform temperature {T x ) nonscattering medium en- 
closed in a black cold cylinder with unity radius and height 
equivalent to one diameter. The net surface heat flux q s (or, 
in this particular case, heat flux q, at r = 1), by three different 
solutions is shown in the figure. The YIX results are very 
close to the “exact” solution by Dua and Cheng. 1 '' 1 The max- 
imum differences between YIX results and the exact solution 
are 3.2, 1.9, and 9.6% for r = 5.0, 1.0, and 0.1, respectively. 
The P3 solution by Sun, 2 " which is similar to that of Menguc 
and Viskanta, 15 is reasonably accurate for r = 1, except at 
large z, where it overpredicts. At low optical thickness, P3 
has a large deviation from the exact solution. The finite vol- 
ume method 4 results are nearly the same as the exact solutions 
and are not shown in the figure. The calculation time for the 
YIX method increases as the optical thickness increases. The 
CPU times for r = 5.0. 1.0, and 0.1 are 335, 195, and 32 min 
respectively, with 0.001 as the first integration point and S16 
as the angular quadrature. 

Experimental cylindrical furnace data were available for 
comparison with several numerical solutions. The 0.9-m diam 



Fig. 2 NondimensionaJ radiant heat flux on the wall of cylindrical 
enclosure containing nonscattering medium at three different optical 
thicknesses. 



AXIAL DISTANCE, m 


Fig. 3 Comparison of local surface heat fluxes at the side wall of a 
cylindrical furnace with experimental data. 

and 5-m length water-cooled Delft furnace data were obtained 
from Jamaluddin and Smith 16 and Wu and Fricker. 21 The 
radiative medium for this experiment was modeled as gray, 
nonscattering with constant extinction coefficient of 0.3 m' 1 . 
The boundaries are at 425 K with an emissivity of 0.8, except 
at the inlet and exit plenums, where the boundaries are treated 
as black surfaces at 300 K. Figure 3 shows the measured wall 
heat flux and calculations by various methods. All the nu- 
merical solutions correctly predict the location of maximum 
heat flux. However, the P3 method 15 seriously underpredicts 
the value of maximum heat flux, while other methods are 
reasonably close to each other. In the legend in Fig. 3, YIX/ 
Sn represents the YIX quadrature using 0.001 as the first 
integration point and an angular quadrature using the Sn dis- 
crete ordinates set. It is interesting to note that the YIX/S4 
and finite volume 4 methods produce nearly identical results. 
The higher-order quadrature (YIX/S16) predicts a higher peak 
flux than all other methods at an axial distance 1 m from the 
furnace base. Compared with results from the YIX and finite 
volume methods, the S4 result 22 is lower for almost the entire 
length of the furnace. It is acknowledged 4 ' 5 * 22 that the mea- 
sured values cannot be reproduced even with more accurate 
numerical methods. One of the reasons for the discrepancy 
is believed to be the assumption of a uniform extinction coef- 
ficient. This shows the necessity for treating the medium as 
nonhomogeneous in radiation transfer calculation for com- 
bustion systems. 

For the experimental furnace problem, the YIX run times 
and maximum errors for using different integration points and 
angular quadrature orders are shown in Table 1. Note that 
the SI 6/0.001 case is used as the baseline for error comparison. 
As expected, with the higher-order discrete ordinates set, the 
CPU time increases in proportion to the corresponding N w . 
While the smaller first integration point used in the YIX quad- 
rature reduces the integration error, it also increases the CPU 
time. 

As part of the effort to model a turbulent jet diffusion 
sooting flame, the radiation heat transfer within such a system 
is treated rigorously. Flame temperature (Fig. 4) and soot 
volume fraction data are obtained from numerical results for 
a turbulent ethylene jet diffusion flame in quiescent air with 
a nozzle diameter of D = 0.58 mm and a fuel flow rate of 
3.96 cm 3 /s. 2J Since the extinction coefficient of soot aggregates 
can be approximated by that of Rayleigh spheres, 24 it can be 
calculated from volume fraction with known soot refractive 
index and wavelength. We chose m = 1.7 - /0.7 and 0.5- 
ptm wavelength for this study. We are not claiming that the 
medium should be assumed gray. Rather, we do this so that 
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Table 1 Comparison of CPU time and maximum error for the example problem 

in Fig. 3 




S 16 


S4 

The first YIX 

CPU 

Maximum 

CPU 

Maximum 

integration point 

time, s 

error. % 

time, s 

error, % 

0.001 

270 

0 (base case) 

24.8 

4.2 

0.01 

37.4 

4.7 

3.8 

5.9 


A 3 x 17 grid system is used in all calculations on a Sun workstation. 



D 


Fig. 4 Temperature (K) distribution of a turbulent jet diffusion flame. 



D 


Fig. 5 Extinction coefficient (m~‘) distribution of a turbulent jet 
diffusion flame. 

we can focus on effects of nonhomogeneity and scattering and 
be more computationally efficient. The local extinction coef- 
ficient, shown in Fig. 5, varies from as low as 10" 13 m ~ 1 to 
about 10m' 1 within the whole computational domain, which 
is much larger than the region covered in the figures. The 
computational domain is between 0 < rID ^ 130 and 0 ^ 
zID ^ 500 using a 65 x 50 grid. The boundaries are assumed 
to be black at 300 K. A high-order integration quadrature is 
used: the first YIX integration point of 0.001 and the S16 
discrete ordinates set. 

We first calculate the nonscattering case. Based on the data 
of Figs. 4 and 5, the calculated normalized radiative heat flux 



Fig. 6 Normalized radiative flux divergence contour and radial heat 
flux (at r/D = 31) for turbulent diffusion flame with nonhomogeneous 
extinction coefficient and scattering albedo equals 0. 



Fig. 7 Normalized radiative flux divergence contour and radial heat 
flux (at rID = 31) for turbulent diffusion flame with homogeneous 
extinction coefficient and scattering albedo equals 0. 


divergence distribution is plotted along the r-z axis (the left 
part of Fig. 6). The peak value (=370) of the normalized flux 
divergence occurs near the flame center, i.e., at rID ~ 0 and 
z!D = 155. The radiative heat flux is normalized with respect 
to the blackbody emissive power at 1000 K. In Fig. 7 the same 
temperature data are used, but the extinction coefficient is 
treated as a constant. An equivalent mean extinction coeffi- 
cient is obtained by averaging the local extinction coefficients 
of all the volume elements whose temperatures are higher 
than 300 K. As shown in Fig. 7, the peak flux divergence is 
still at the flame center, but its value is reduced to about 49 
from 370, by a factor of more than 7. In the right parts of 
Figs. 6-10, the corresponding q r at r!D = 31 for different 
conditions is also plotted. The homogeneous case (Fig. 7) has 
smooth q r distribution due to its constant extinction coeffi- 
cient. On the other hand, nonhomogeneous cases have a sharp 
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Fig. 8 Normalized radiative flux divergence contour and radial heat 
flux (at rlD = 31) for turbulent diffusion flame with non homogeneous 
extinction coefficient and scattering albedo equals 0.2. 



Fig. 9 Normalized radiative flux divergence contour and radial heat 
flux (at HD = 31) for turbulent diffusion flame with nonhomogeneous 
extinction coefficient and scattering albedo equals 0,5. 



r__ Jl_ 

D oT 

Fig. 10 Normalized radiative flux divergence contour and radial heat 
flux (at r!D - 31) for turbulent diffusion flame with nonhomogeneous 
extinction coefficient and scattering albedo equals 0.8. 


peak q, . It is evident that the effect of neglecting radiative 
property nonhomogeneity is very significant. To predict the 
flame behavior accurately with strong radiation effect and 
nonuniform soot distribution, it is necessary to model the 
radiation heat transfer by considering the local radiative prop- 
erty variation. 

Figures 8-10 show the results under the same conditions 
as Figs. 4 and 5, but with isotropic scattering added where 
the scattering albedo equals 0.2, 0.5, and 0.8, respectively. 
Recent analysis 24 and measurements 25 both show that the scat- 
tering effect can be very significant, since soot particles are 
agglomerated in flames. Figure 8 shows a nearly identical flux 
divergence distribution as a nonscattering case (Fig. 6), except 
that the peak value is reduced to about 330 from 370 and the 
radial heat flux also decreases. As the scattering albedo in- 
creases, both radiative flux and its divergence decrease. Note 
that the scattering albedo may be a function of the position 
in the flame, since at later stages of the combustion (or at a 
higher position above the jet nozzle), the soot aggregates may 
grow larger, which in turn increases the scattering albedo. In 
the current numerical scheme using the YIX method, it is 
very easy to incorporate the distribution of the scattering 
albedo. However, since this data is not yet available, the 
calculation will be left for future study. 

Conclusions 

Numerical solutions of radiative heat transfer within a finite 
axisymmetric cylindrical enclosure involving homogeneous and 
nonhomogeneous media are treated rigorously. The use of 
the YIX quadrature and discrete ordinates is Shown to be able 
to solve complicated multidimensional homogeneous as well 
as nonhomogeneous problems accurately. Numerical results 
for soot radiation from a turbulent jet diffusion flame show 
that neglecting the nonhomogeneity in radiative property can 
cause significant errors. They also show that the magnitude 
of scattering albedo has a significant effect on the radiative 
heat flux. It should be pointed out that the current numerical 
solution can also easily include nonuniform distribution of 
any other radiative properties, such as the scattering albedo. 
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